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Summary 

This report studies the effects of fractional dynamics in 
chaotic systems. In particular, Chua’s system is modified to 
include fractional order elements. Varying the total system 
order incrementally from 2.6 to 3.7 demonstrates that systems 
of “order* * less than three can exhibit chaos as well as other 
nonlinear behavior. This effectively forces a clarification of the 
definition of order which can no longer be considered only by 
the total number of differentiations or by the highest power of 
the Laplace variable. 

Introduction 

It is well known that chaos cannot occur in continuous-time 
systems of order less than three. This assertion is based on the 
usual concepts of order, such as the number of states in a 
system, the highest power of the Laplace variable, s , in the 
system, or the total number of separate differentiations or 
integrations in a system. Unfortunately, these concepts of order 
do not directly relate to systems having fractional order com- 
ponents. The purpose of this report is to demonstrate that 
systems of order less than three, as defined in the usual way, can 
still display chaotic behavior. The next section provides a brief 
review of fractional calculus. Useful approximations for these 
fractional operators follow. Finally, an example is given which 
demonstrates that systems of order less than three can display 
chaos. This is both shown experimentally via simulations and 
predicted analytically using the describing function method. 

Review of Fractional Operators 

The idea of fractional integrals and derivatives has been 
known since the development of the regular calculus, with the 


first reference probably being associated with Leibniz in 1695 
(Oldham and Spanier, 1974, page 3). Although not well known 
to most engineers, the fractional calculus has been considered 
by prominent mathematicians (Courant and Hilbert, 1953) as 
well as the “engineers” of the operational calculus (Heaviside, 
1971; and Bush, 1929). In fact several textbooks written before 
1 960 have some small section on fractional calculus (Goldman, 
1949; Holbrook, 1966; Starkey, 1954; Carslaw and Jeager, 
1948; Scott, 1955; and Mikusinski, 1959). An outstanding 
historical survey can be found in Oldham and Spanier (1974) 
who also give what is unquestionably the most readable and 
complete mathematical presentation of the fractional calculus. 
Other bound discussions of the area are given by Ross (1975), 
McBride (1979), and McBride and Roach (1985). Unfortu- 
nately, many of the results in the fractional calculus are given 
in the language of advanced analysis and are not readily 
accessible to the general engineering community. 

Many systems are known to display fractional order dynam- 
ics. Probably the first physical system to be widely recognized 
as one demonstrating fractional behavior is the semi-infinite 
lossy (RC) line. The current into the line is equal to the half- 
derivative of the applied voltage; that is, the impedence is 

V(s) = -L/(i) 

Vi 

Although this system was studied by many, Heaviside (1971) 
considered it extensively using the operational calculus. He 
states “there is a universe of mathematics lying in between the 
complete differentiations and integrations” and that “fractional 
(operators) push themselves forward sometimes, and are just as 
real as the others.” Another equivalent system is the diffusion 
of heat into a semi-infinite solid. Here the temperature looking 
in from the boundary is equal to the half integral of the heat rate 
there. Other systems 'that are known to display fractional order 
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dynamics are viscoelastic systems (Bagley and Calico, 1991, 
Koeller, 1984;Koeller, 1986; Skaar, Michel, and Miller, 1988; 
Lopez-Marcos, 1990 ); colored noise (Mandelbrot, 1967); elec- 
trode-electrolyte polarization (Ichise, Nagayanagi, andKojima, 
1971; Sun, Onaral, and Tsao, 1984); dielectric polarization 
(Sun, Abdelwahab, and Onaral, 1984); boundary layer effects 
in ducts (Sugimoto, 1991); and electromagnetic waves 
(Heaviside, 1971). Because many of these systems depend 
upon specific material and chemical properties, it is expected 
that a wide range of fractional order behaviors are possible 
using different materials. 

Two commonly used definitions for the general fractional 
differintegral are the Grunwald definition and the Riemann- 
Liouville definition (Oldham and Spanier, 1974). The Ri- 
emann-Liouville definition of the fractional integral is given 
here as 


d q f 

dt q 


— 1' 
n-<?) Jo 


m 

(f-T )* +1 




q < 0 


where q can have noninteger values, and thus the name frac- 
tional differintegral. Notice that the definition is based on 
integration and more importantly is a convolution integral for 
q < 0. When q > 0, then the usual integer nth derivative must be 
taken of the fractional (q - n)th integral, and yields the frac- 
tional derivative of order q as 


d q f _d*_ d q ~ n f 
dt q ” dt n dt q ~ n 


q > 0 and n an integer > q 


This appears so vastly different from the usual intuitive defini- 
tion of derivative and integral that the reader must abandon the 
familiar concepts of slope and area and attempt to get some new 
insight (which still remains elusive). This is discussed further 
in Lorenzo, C.F.; and Hartley, T.T.: On Conceptualization, 
Initialization, and Applications in Fractional Calculus (to be 
published). 

Fortunately, the basic engineering tool for analyzing linear 
systems, the Laplace transform, is still applicable and works as 
one would expect; that is, 


l\^^\= s q L {m} 


Amazingly enough, one of the most difficult obstacles in the 
practical application of the fractional calculus is the initial 
condition problem. As long as a given system is at rest, at the 
zero equilibrium, at time zero, the fractional initial value 
problem is readily solved using standard Laplace transform 
methods (all initial condition terms are zero). Unfortunately, 
the fractional derivative operator starts rather abruptly at time 
zero, so that any nonzero initial value for a function will appear 
as a discontinuity and translate directly into a t r term, which 
has an annoying singularity at time zero using the appropriate 
power, r. This is not necessarily a problem, unless the desired 
initial value is not infinity. Bagley (1988) addresses this prob- 
lem by creating a modified fractional derivative operator that 
essentially subtracts out the singularity. The problem is further 
studied by Hartley, T.T.; and Lorenzo, C.F.: Insights Into the 
Fractional Initiative Value Problems (to be published) by 
relating it back to the semi-infinite line problem. 

Bagley (1988) has also extended the initial value problem to 
fractional state space systems. Here the idea of state no longer 
gives all past and future knowledge of the system behavior via 
some stored pseudo-energy. In fact, the number of these frac- 
tional states is somewhat arbitrary and dependent only upon 
what the user has chosen as the base fractional derivative. 

Understanding the possible dynamic behavior of linear frac- 
tional order systems is fundamental to the development of 
future applications. Progress in this area has been fairly slow, 
however, since there was no known general fractional order 
impulse response with which to perform convolution. 
Recently, Bagley (1988) has shown that the impulse responses 
of fractional order systems are related to the Mittag-Leffier 
function (Erdelyi, et al. 1955), which is effectively the frac- 
tional order analog of the exponential function. With this 
knowledge, it has been possible to better clarify the time 
responses associated with fractional order systems. Impulse 
responses, step responses, and initial condition responses for 
some general fractional order systems can be found in Hartley, 
T.T.; and Lorenzo, C.F.: The Solution to a General Linear 
Fractional Order Initial Value Problem (to be published). 


[ dt q J *=o 

where n is an integer such that n - 1 < q < n (Oldham and 
Spanier, 1974). If the initial conditions are considered to be 
zero, this formula reduces to the more expected and comforting 
form 


The Concept of System Order 

As the concept of “order” is central to the understanding of 
fractional systems, some discussion of this concept now fol- 
lows. In this discussion, it will be assumed that the systems 
being considered are single-input-single-output, that their 
representations are minimal in the usual sense (Kailath, 1980), 
and that they are linear. 


dr l ~ k m 

dt q -'~ k 


, for all q 


Jf=o 
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Mathematical order is defined as the highest derivative 
occurring in a given differential equation. The concept of 
mathematical order is applicable to both ordinary and fractional 
differential equations. Normally, when the word “order” is 
used without a qualifier, it implies the meaning of mathematical 
order. 

For linear dynamic systems that are described by ordinary 
differential equations (i.e., of integer mathematical order), the 
system mathematical order implies, or is equivalent to, the 
following: 

(1) The highest derivative in the ordinary differential 

equation 

(2) The highest power of the Laplace variable, 5 , in the 

characteristic equation 

(3) The number of initializing constants required for the 

differential equation 

(4) The number of singularities in the characteristic 

equation 

(5) The length of the state vector 

(6) The number of energy storage elements 

(7) The number of independent spatial directions in which 

a trajectory can move 

(8) The number of devices that add 90° sinusoidal steady 

state phase lag 

(9) The number of devices that retain some memory of the 

past 

The utility of the definition of mathematical order is that it 
infers all the system characteristics for systems with only 
integer order components. 

Thus the benefit of having a definition for order for linear 
ordinary differential equations is that it allows a direct under- 
standing of the behavior of a given dynamic system. Unfortu- 
nately, for fractional differential equations, the order of the 
highest derivative does not infer (or is not equal to) all of the 
previously mentioned properties. Indeed, the most important 
characteristic of order in integer order ordinary differential 
equations is probably item (3) in the previous list (i.e., it 
dictates the number of initializing constants which together 
with the differential equations allow prediction of the future 
behavior). In systems terminology, this information provides 
the initial “state” of the system being analyzed. Clearly, the 
order of the highest derivative in a fractional differential equa- 
tion does not have this property, nor does it predict the associ- 
ated number of energy /memory elements associated with the 
fractional differential equation, nor does it infer the number of 
integrations (even fractional) required to solve or simulate the 
given fractional differential equation. Thus the issue of order 
and the information required together with the fractional differ- 
ential equation to predict future behavior is fundamental and is 
expected to be treated in detail at a later time. 


Approximation of Fractional Operators 

The standard definitions of the fractional differintegral do 
not allow direct implementation of the operator in time-domain 
simulations of complicated systems with fractional elements. 
Thus, in order to effectively analyze such systems, it is neces- 
sary to develop approximations to the fractional operators 
using the standard integer order operators. In the work that 
follows, the approximations are effected in the Laplace s- 
variable. The resulting approximations provide sufficient accu- 
racy for time domain hardware implementations. 

Some work has been done in this area already, but it has riot 
been highly organized. Oldham and Spanier (1974) and Piche 
(1992) give several discrete-time approximations based on 
numerical quadrature. In continuous time, engineers have used 
network theory approximations (Carlson and Halijak, 1964; 
Steiglitz, 1964; Carlson and Halijak, 1961 ; and Halijak, 1 964). 
More recently Oldham and Spanier ( 1974), Ichise, Nagayanagi, 
and Kojima, 1971; and Charef, et al. (1992) have developed 
other network theory approximations. Even more recently, a 
discrete-time fractional calculus has been developed similar to 
the theory of linear multi step methods for numerical integration 
(Lubich, 1985, 1986, 1988a, and 1988b). 

The approximation approach taken here is that of Charef, et 
al. (1992). Basically the idea is to approximate the system 
behavior in the frequency domain. This is done for a given q by 
creating an approximation with Bode magnitude response roll 
off of 20 times q db/dec, which will consequently have a phase 
shift of approximately 90 times q degrees over the required 
frequency band. This approximation is created by choosing an 
initial breakpoint (the low frequency accuracy limit of the 
approximation), the allowable error in db 1 s, and the number of 
r-plane poles in the approximation. The high frequency limit of 
the usable bandwidth can be varied by changing the allowable 
error and the number of poles. Thus an approximation of any 
desired accuracy over any frequency band can be achieved. 
Table I gives approximations for \!s q with q = 0.1 to 0.9 in 
steps of 0.1. These were obtained by trial and error and are 
reasonably good from 0.01 to 100 rad/sec. These approxima- 
tions are used in the study that follows. 


A Fractional Chua System 

Chua’s system is well known and has been extensively 
studied. The particular form to be considered here was pre- 
sented by Hartley (1989) and used further for the study of 
Hartley and Mossayebi (1993). This system is different from 
the usual Chua system in that the piecewise-linear nonlinearity 
is replaced by an appropriate cubic nonlinearity which yields 
veiy similar behavior. It is represented in state space form as 
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y = x - y + Z 



It is studied here in two different system representa- 
tions as discussed in the following sections. In all cases studied, 
P is defined to be 100/7 and a is allowed to vary. 

State Space Configuration 

To study the effect of fractional derivatives on the dynamics 
of this system, the state space configuration (fig. 1(a)) was 
considered first. Here, the vector derivative was replaced by a 
vector fractional derivative as follows: 


prohibited a timely calculation of any exponents but the first. 
Since the order of this system was greater than three, these 
calculations were not pursued. In all cases, the one positive 
exponent clearly indicated that the system was behaving cha- 
otically. The numerical simulations also indicated that the 
lower limit of the vector fractional derivative q was between 
0.8 and 0.9 for this system to remain capable of generating 
chaos. The lowest value obtained for mathematical order to 
yield chaos was 2.7 using the q = 0.9 fractional vector deriva- 
tive. No upper limit was obtained. Phase plane plots for these 
systems are given in figure 2. 

Feedback Configuration 

The feedback configuration is now considered. To change 
the total system mathematical order, the separated 1 Is in figure 
1 (b) was allowed to change powers, that is. 


1 1 






s l00 y--8y 
dt q i 


Simulations were then performed using q - 0.8, 0.9, 1.0, and 
1.1. The approximations from table I were used for the simula- 
tions of the appropriate <?th integrals. When q < 1, then the 
approximations were used directly. It should further be noted 
that approximations used in the simulations for 1/s*, when 
q > \ 9 are obtained by using 1/s times the approximation for 
1/s*" 1 from table I. 

The results from this state space study verified that chaos 
could indeed occur in a system of mathematical order less than 
3. This was determined by computing the Lyapunov exponents 
for each of the simulations with q = 0.9, 1 .0, and 1.1, using the 
method of Benettin, et al. (1990). Chaos is indicated when any 
of the Lyapunov exponents is greater than zero. These results 
are given in table II where the largest Lyapunov exponents are 
given as a function of system order. In each case, the second 
exponent was near zero. The 2.7 order system approximation 
had an additional six negative exponents which were not listed. 
Also the 3.3 order system approximation was so large that it 


A variety of simulations were performed on the resulting 
systems as discussed subsequently. Here, the approximations 
from table I were used to represent the fractional integral where 
again the approximations for 1/s*, when q > 1 , were obtained 
by using 1/s times the approximation for 1/s*" . 

Bifurcation diagrams for several of these systems are given 
in figure 3. Here, a particular value of q was chosen, and the 
parameter was varied to obtain the particular bifurcation plot. 
These diagrams were generated by simulation using Euler’s 
method and a simulation timestep of 0.001. These were veri- 
fied by further reducing the timestep by an order of magnitude 
with little change in the overall bifurcation structure. To obtain 
these diagrams, the values of the output x- variable were plotted 
whenever its slope changed sign. Although it is believed that 
the bifurcation diagrams are reasonably accurate and are 
sufficiently accurate for this particular study, more correct 
diagrams could possibly be obtained by using more accurate 
approximations of the fractional derivative than those given in 
table I or a more accurate simulation. Observation of the 
bifurcation diagrams indicates behavior similar to that from 
the state space study. For the feedback configuration, decreas- 
ing the power of s shifts the bifurcation diagram to the right as 
a function of a, while the converse is also true. The limits on 
the system mathematical order to have a chaotic response as 
measured from the bifurcation diagrams are approximately 2.5 
< n < 3.8. The overall behavior from the simulation studies is 
summarized in figure 4. 

An advantage to the feedback configuration is that it allows 
easy system analysis using describing functions, as discussed 
in Hartley and Mossayebi (1993). Here the idea is that the 
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frequency response of the linear block in the feedback configu- 
ration is plotted in the Nyquist plane, together with minus one 
over the appropriate describing function of the nonlinearity, as 
in figure 5. The fractional order integral in the loop is handled 
directly by taking the frequency response on the primary 
Riemann sheet and essentially poses no complication or confu- 
sion in application of the describing function approach. In other 
words, the fact that fractional powers of s are present does not 
require any frequency domain approximation as in the time- 
domain simulation; rather the fractional powers of s can be used 
as is in computing the frequency response of the linear block. 
In Hartley and Mossayebi ( 1 993), it is shown that the important 
points from the nonlinearity of this system in the Nyquist plane 
are 

(1) ^e[//(yco))] > -3.5, <9wf//(/co)] = 0, which indicates two 
stable points at x = ±^lo^5. 

(2) ^e[//(/co))] < -3.5, 5w[H(/co)] = which indicates a 
Hopf bifurcation of the stable points of item (1) into a limit 
cycle. 

(3) <Be[//(/a)))] ^ -7, 5w[//(/co)] = 0, which indicates that 
period doubling of the limit cycle of item (2) occurs (this 
progresses into spiral chaos). 

(4) <$e[//(/co))] < -14, 5w[//(/co)] = 0, which indicates merg- 
ing of the spiral chaos into the double scroll behavior. 

Extinction of the double scroll (meaning its disappearance) is 
not directly predicted using the describing function approach, 
but a reasonable approximate value is <£e[//(/co))] < -23, 
5w[//(/co)] = 0. A diagram indicating the usage of the describing 
function is given in figure 5. 

Using these results and varying the power of the integrator in 
the loop allowed a theoretical prediction of the simulation re- 
sults of figure 4. These theoretical results are given in figure 6. 
It should be noted that the qualitative features are very well 
predicted using the describing function approach, and that the 
quantitative results are reasonably close. Furthermore, for 
mathematical system order less than approximately 2.85, the 
describing function approach predicted the appearance of a 
stable and unstable limit cycle as a increased (via an apparent 
saddle node bifurcation). These limit cycles coexist with each 
of the stable fixed points. Eventually, as a increased further, the 
unstable cycles merged with the stable fixed points via a 
subcritical Hopf bifurcation, leaving unstable fixed points. 
This entire process basically became a supercritical Hopf 
bifurcation for mathematical order greater than 2.85. This was 
then verified in the simulations with this bifurcation structure 
occurring for mathematical system order less than approxi- 
mately 2.75. In fact, for the mathematical order equal to 2.6, the 
simulation showed the points at * = ±Vo 5 to be stable and 
each coexisting with spiral chaos. It is a true testament to the 
utility of the describing function approach that it could predict 
the behavior of this system as accurately as it does. 


Concluding Remarks 

This report has introduced the idea of fractional derivatives 
from the dynamic systems viewpoint. It has been demonstrated 
that the usual idea of system order must be modified when 
fractional derivatives are present. The usual approach of calcu- 
lating the mathematical system order by determining the 
highest derivative in the system does not work in this situation. 

It has been further demonstrated that chaos, as well as the 
other usual nonlinear dynamic phenomena, can occur in sys- 
tems with mathematical order less than three via Chua’s 
system. This is surprising given the usual nonlinear system 
paradigms concerning chaos and order. It is not clear at this 
point whether the chaos in fractional order systems should be 
characterized differently than chaos in regular integer order 
systems. 

It should be noted that the describing function approach 
usually requires at least -1 80° of phase shift in the linear part 
of the feedback loop to ever predict Hopf bifurcations, and 
consequently chaos, for memoryless nonlinearities. Because 
the linear part can be a nonminimum phase transfer function, 
it is further conjectured that chaos can occur in systems with 
mathematical order less than three and probably less than one. 
Furthermore, the feedback configuration indicates that, as long 
as the linear part of the loop has at least -180° of phase shift, 
the possibility of chaos in the system depends primarily on the 
nonlinearity and how its particular describing function be- 
haves. 

As has been demonstrated, the idea of fractional derivatives 
requires one to reconsider dynamic system concepts that are 
often taken for granted. Some of these concepts have been 
discussed in this report. Some others that require much further 
consideration are the concept of Lyapunov exponents for 
fractional states, the use of fractional states in which to embed 
attractors, and the relationship between fractional order and 
fractal dimension. 

National Aeronautics and Space Administration 

Lewis Research Center 

Cleveland, Ohio 44135, June 30, 1995 
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TABLE I.— FRACTIONAL OPERATORS WITH APPROXI- 
MATELY 2 db ERROR FROM w = HU 2 TO 10 2 rad/sec 


1 220.4/ + 5004/ + 5038.5 2 + 234.5s + 0.4840 

0.1 

s 

/ + 359.8/ + 5742/ + 4247./ + 147.75 + 0.2099 

1 

60-95/ + 816.9s 3 + 582.8s 2 + 23.24s + 0.04934 

0.2 

5 

/ + 134.0/ + 956.5j 3 + 383.5/ + 8.9535 + 0.01821 

1 

23.765 4 + 224.95 3 + 129.\s 2 + 4.7335 + 0.01052 

0.3 

5 

/ +64.515 4 + 2S2.25 3 + 63.615 2 +1.1045 + 0.002267 

1 

25.00s 4 + 558.5s 3 + 664.2s 2 + 44.15s + 0.1562 

0.4 

5 

/ + 125.6 5 4 + 840.65 3 + 317.25 2 + 7.4285 + 0.02343 

1 

15.97s 4 + 593.2s 3 + 1080s 2 + 135.4s + 1 

0.5 

5 

/ + 134.35 4 + 10725 3 + 543. 45 2 + 20.105 + 0.1259 

1 

8.5795 4 + 255.6 5 3 + 405.3 s 2 + 35.93 5 + 0.1696 

0.6 

5 

5 5 + 94.225 4 + 472.95 3 + 134.85 2 + 2.6395 + 0.009882 

1 

5 .4065 4 + 177.65 3 + 209.65 2 + 9.1975 + 0.01450 

0.7 

5 

/ + 88.12s 4 + 279.2s 3 + 33.30s 2 + 1.927s + 0.0002276 

1 

5.235s 3 + 1453s 2 + 5306s + 254.9 

0.8 

5 

5 4 + 658. l5 3 + 57005 2 + 658.2 5 + 1 

1 

1.7665 2 +38.275 + 4.914 

0.9 

5 

5 3 + 36.155 2 + 7.7895 + 0.01000 


TABLE II — LARGEST LYAPUNOV EXPONENTS FOUND IN THE 
STATE SPACE CONFIGURATION FOR q = 0.9, 1.0, AND 1.1 
WHICH GIVES A TOTAL SYSTEM MATHEMATICAL 
ORDER OF 2.7, 3.0, AND 3.3, RESPECTIVELY 


Mathematical 
system order 

Order of system 
approximation 

a-used 

Exponents 


*2 

*3 

2.7 

9 

12.75 

0.338 

-0.000201 

-0.132 

3.0 

3 

9.50 

0.248 

-0.00412 

-3.07 

3.3 

18 

7.00 

0.318 

(a) 

(a) 


a These values were not calculated. 
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Figure 2.— Phase plane projections for the state space configuration of Chua's system, 
(f = 200, AT = 0.05.) (a) Total mathematical system order is 3.0, ot = 9.5. (b) Total 
mathematical system order is 3.0, a = 9.5. 
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Figure 2.— Concluded. Phase plane projections for the state space configuration of Chua's 
system, (f = 200, AT = 0.05.) (i) Total mathematical system order is 3.3, a = 7.0. 
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- 1.0 1 

12.0 


Figure 3. — Bifurcation diagram for the feedback configuration of Chua‘s system; 
maximum and minimum of x plotted against a. (a) Fractional integral of order 0.6 
total mathematical system order 2.6. (b) Fractional integral of order 0.7, total 
mathematical system order 2.7. 






Figure 3.— Continued. Bifurcation diagram for the feedback configuration of Chua's 
system; maximum and minimum of x plotted against a. (c) Fractional integral of order 
0.8, total mathematical system order 2.8. (d) Fractional integral of order 0.9, total 
mathematical system order 2.9. 











Figure 3.— Continued. Bifurcation diagram for the feedback configuration of Chua's 
system; maximum and minimum of x plotted against a. (e) Fractional integral of order 
1 .0, total mathematical system order 3.0. (f) Fractional integral of order 1 .1 , total 
mathematical system order 3.1. 
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Figure 3.— Continued. Bifurcation diagram for the feedback configuration of Chua's 
system; maximum and minimum of x plotted against a. (g) Fractional integral of order 
1 .2, total mathematical system order 3.2. (h) Fractional integral of order 3.1 , total 
mathematical system order 3.3. 
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System mathematical order 


Figure 4. — Bifurcation diagram in the a versus system mathematical order plane based on 
simulation studies of the fractional Chua system. Note that the saddle-node and sub- 
critical Hopf merge at 2.75. 



Figure 5. — Nyquist plane plot showing the frequency response of the 
linear part of figure 1 (set of curved lines) for various q and a = 9.5; 
and the describing function of the nonlinearity (solid line on real- 
axis from 3.5 to -20 shown). 
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